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Abstract 

This work is aimed to develop a new class of methods for the BGK model of the Boltzmann 
equation. This technique allows to get high order of accuracy both in space and time, teoretically 
without CFL stability limitation. It's based on a Lagrangian formulation of the problem: infor- 
mation is stored on a fixed grid in space and velocity, and the equation is integrated along the 
characteristics. The source term is treated implicitly by using a DIRK (Diagonally Implicit Runge 
Kutta) scheme in order to avoid the time step restriction due to stiff relaxation. In particular 
some L-stable schemes are tested by smooth and Riemann problems, both in rarefied and fully 
fluid regimes. Numerical results show good accuracy and efficiency of the method. 

1 Introduction 

The Kinetic Theory of gases is based on the Boltzmann equation (BE) for 
the distribution function f(t, x, v), which depends on the time t, the physical 
space coordinates vector x and the microscopic velocity space coordinates 
vector v,and provides an accurate description of rarefied gas flows. The rar- 
efied regime arises once the Knundsen number Kn = X/L = O(l), where A 
is the mean free path and L is a certain characteristic length. Because of its 
nonlinearity and multi-dimensionality, it is very difficult to find analytical 
solutions of BE. In most pratical cases, numerical methods have to be used. 
The Boltzmann equation has a wide range of applications. It is being used to 
study flows in external atmosphere, in particular flows around spacecrafts. 
This kind of applications are characterized by so large Knundsen numbers 
that the kinetic approach gets highly demanded. Recent applications con- 
cerns flows in microchannels and Micro Electro Mechanical Systems (MEMS) 
|15] , Here the gas works usually at standard conditions, whereas the mean 
free path lies on submicron scale, like the characteristic length of the system. 

The numerical solution of the Boltzmann equation leads us to be involved 
in a very challenging topic. This is due mainly to the high dimensionality of 
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the problem (in the fully 3D case the distribution function depends on seven 
independent variables), the nonlinear collision term, and the requirement to 
preserve the collision invariants at a discrete level. Several strategies have 
been developed to tackle the problem, each one being suitable for appropriate 
circumstances. The numerical methods may be grouped in two main classes. 
One is the class of probabilistic methods, like Direct Simulation Monte Carlo 
(DSMC) method described by Bird 0], or recently by Pareschi and Russo 
[30] . The second one concerns deterministic methods, see Buet [6], Ohwada 
|24] . Rogier and Schneider |36| ; for a review on both classes of method see 
|31j . The most demanding numerically part concerns the computations of 
the collision term. The peculiarity of the deterministic methods lies in the ca- 
pability of making us able to get very accurate solutions about special cases. 
Belonging to the deterministic class, spectrally accurate methods have been 
proposed (see, for example, |32| [TTJ [12]), though the computational cost is 
high. The Monte Carlo methods, on the other hand, allow computations 
of the collision term in very efficient way. They are easy-fitting, and can 
handle physical problems arising from strongly nonequilibrium conditions. 
Due to the statistical noise, accurate solutions can be pursued just by many 
average steps. Therefore, Monte Carlo methods are not efficient for systems 
near equilibrium, whereas deterministic methods are pretty demanded. An 
exaustive example may be the computation of gas flows in MEMS, see |14j . 
particularly characterized by low velocities. Pareschi and Caflisch |25| re- 
cently proposed an alternative approach, by modifing the DMCS. However 
the trouble concerning the statistical noise is not fully worked out yet. 

The computation time can be reduced by considering simplified models 
of the Boltzmann equation, like the BGK model, introduced by Bhatnagar, 
Gross and Krook [3]. BGK model shows good properties. One most of 
all in particular, the convergence to the Euler equation once the Knundsen 
number approaches to zero. In some cases it works well enough also far 
from the equilibrium, see |16j . This model has been extensively theoretically 
investigated |26| . and many numerical computations have been carried out 
in order to validate its properties (|10|.|40|). Some interesting applications 
of BGK model are described in [1] and [2]. Moreover BGK model is also 
applied to particular flows in nanostructure |17] , The classical schemes tipi- 
cally suffer the inefficiency due to the stiffness arising from the relaxation 
time getting smaller and smaller. Standard schemes require the solution of 
nonlinear systems derived by discrete formulations of the problem in implicit 
form. Recently this problem has been circumvented by Puppo and Pieraccini 
|27] by an implicit formulation of the relaxation term that can be explicitly 
computed. In this work the problem is formulated in semilagrangian form. 
The discrete form is worked out on a fixed square grid into the phase space, 
whereas the time integration is performed along the characteristics. The 
time integration is implemented by treating implicitly the relaxation term 
with a technique similar to the one used in |27j . However, because of the La- 
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grangian nature, the usual CFL stability restriction does not apply. In order 
to perform the time evolution along the characteristics it is needed to recon- 
struct all components of the distribution function for each space grid node. 
This task is solved by using high order pointwise WENO reconstruction de- 
scribed in [7]. This class of methods join together high order accuracy and 
high efficiency because the CFL limitation is definitely avoided. Moreover, 
they work up to the fluid dynamic regime, though they suffer from accuracy 
loss when the Knudsen number gets small. Although the schemes are not 
strictly conservative, numerical tests show that the conservation errors are 
very small for smooth flows. A fully conservative version of the schemes can 
be constructed. This paper is organized as follow: after this introduction, in 
section 2, we recall the BGK models and some of its properties. Section 3 is 
devoted to the detailed description of the not conservative schemes, whereas 
in Section 4 we present the result of some numerical tests. In Section 5 
we present the conservative version of the method, and show some results 
pursued by it. Finally, in the last section we give some conclusions. 



2 The BGK model 

The BGK model, introduced by Bhatnagar, Gross and Krook [3], is a sim- 
plification of the Boltzmann equation where the collisions are modeled by a 
relaxation of the distribution function f(t, x, v) towards the Maxwellian. It 
consists of the following initial values problem 

9fit £ ,V) + « • V,/(f , x, v) = l -(M [f(t, x, v)] - f(t, x, v)), 

/(0,x,v) = / (x,v) t>0, xeR d ,veR*. 

Here d > 1 and N > d denote the dimensions of the physical and ve- 
locity spaces respectively. M[f] = M(v;{p,u,T}) is the local Maxwellian 
computed by the moments of the distribution function /(i,x, v) 

M ^te"- r » = (^pv72-p(- ! W ! )- < 2 ' 2 > 

where p = p(x,t), u = u(x, t) and T = T(x, i) denote the macroscopic 
fields, namely: density, mean velocity and temperature, which are related to 
the function / as follows. Let </>(v) = (1, v, l/2v 2 ) T denote the vector of the 
collision invariants of the distribution function /(i,x, v). The moments are 
given by 

(p,pu,Ef = (f^(v)), (2.3) 

where 

(g) = [ giy) dv, g ■. r n m- r. (2.4) 
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The quantity E(x, t) is the total energy and it is related to the temperature 
by the internal energy e(x, t) 

e(x,t) = -RT(x,t), pe = E--pu 2 . 
2.1 Conservation and Entropy Principle 

Conservation laws for the macroscopic fields are regained by (|2.1|) upon mul- 
tiplication by <j) and integration in velocity: 



dt 

d(fv 



W + V x ..(v / )=0, 



dt 



+ V* • (v ® v/ > = 0, (2.5) 



^ + V,<ivVH0. 

The Entropy Principle, sometimes called H-theorem, holds also for the BGK 
models, like for the Boltzmann equation |8] 

9(/ g° g/) +V x -(v/log/)<0, V/(t,x,v)>0. (2.6) 

Once the equilibrium has been established the equality holds, since the dis- 
tribution function is the Maxwellian. 

BGK models are generally implemented by using N = 3, which means that 
the system is a monoatomic gas with three translational degrees of freedom. 
When d = 1 the computation time can be reduced by using the approach 
proposed in pQ, preserving the properties of the gas. In such case, one uses 
the induced cylindrical symmetry in the velocity space, making the problem 
properly one dimensional in space and two dimensional in velocity. As in 
|10| . in this work we choose d = 1, N = 1, since our task is to present the 
methods and testing their properties in simple cases. Thus, the system is a 
gas with one degree of freedom and the integral in velocity space are evalu- 
ated in M. 

The relaxation time of the BGK model is defined by 

r- 1 = cpT 1 ^; (2.7) 

in this definition 5 is the exponent of the viscosity law of the gas (see [9]). 
The constant c is defined by c = RT^ e j/p re j, where p re j is the viscosity at 
the reference temperature T Te f. 

In pQ the authors use expression (|2.7p with 5 = 0, t~ 1 = Ap. 

For clarity of exposition, in this work we assume that the collision fre- 
quency is constant. We shall indicate in a remark how to incorporate in the 
schemes the dependence of the relaxation time on the conservative moments. 
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We reformulate the BGK equation in non-dimensional form, in such a way 
the relaxation time plays the role of the Knudsen number. 

In the BGK model, with a single relaxation time, the transport coeffi- 
cients in the fluid regime (Kn <C 1) are not correctly predicted. For example, 
the Prandtl number Pr is equal to 1, whereas the correct value is |. Be- 
cause of this limitation, several BGK-like models have been proposed (see, 
for example, Struchtrup |39j . Shakov |35j . Liu )22| . Hollway [TO], Bouchut 
and Perthame [5]). For simplicity we will introduce our class of numerical 
methods using the classical BGK model, remarking that the application to 
more sofisticated models is possible as well. 



3 Description of the method 
3.1 A basic first order scheme 

The numerical scheme for the solution of Eq. (|2.ip is based on the charac- 
teristic formulation of the problem (|2.ip . 



df(t,X,v) 1 / , /nfU m f(, \ \ 

= -(M [f(t,x,v)\ - f(t,x,v)), 

at t 

dx (3.1) 
x(0) = x, f(0,x,v) = fo(x,v) t>0, x,v&R. 



Here x becomes a time dependent variable and its equation in (|3.ip can be 
integrated immediately. Hence the BGK model may be presented as follow 

^^ = ±(M[f(t,x,v)]-f(t,x,v)), 

x(t) = x + vt, ^ 3 ' 2 ) 
/(0, x, v) = f (x, v) £ > 0, x,v£R, 

and x is given explicitly. For simplicity we assume constant time step At 
and uniform grid in physical space and velocity domain mesh spacing, Ax 
and Av, respectively and denote the grid points by t n = nAt, X\ = iAx, 
i = 1, . . . , N x , Vj = jAv, j = —N v , . . . , N v , where N x and 2^ + 1 are the 
numbers of grid nodes in space and velocity, respectively. Let ff- denote 
the approximate solution of the problem (|3.2p at time t n in each spatial and 
velocity node. A first order explicit scheme is given by 

Xi = Xi + VjAt, i = l,...,N x , j = —N v , . . . , N v , 
where ( 3 ' 3 ) 
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The numerical solution at time (n 



by m in eq. (J3T3 



, l)At requires f n (xi,Vj), denoted 
The distribution function's value / t " is computed by 



t n + At 
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Figure 1: Characteristics diagram for positive velocity grid node. 

linear interpolation, using the values of the distribution function on the left 
and right nodes of the interval containing Xi = %i — VjAt (see Figure [1]). 

The distribution function for fixed velocity node Vj evolves along the 
characteristics (thick line). The Maxwellian M^[f™] is computed as follow 



M^[/f] 



(2irRTi) N / 2 



exp 



uA 2 



2RT, 



(3.4) 



This formula requires the computation of the discrete moments of ff-. 
This can be done by using a numerical approximation of the integrals com- 
puted in (|2.4p . Following the notation in [23], the discrete velocity grid may 
be denoted by V, which is composed of 2N V + 1 nodes, and (g) can be ap- 
proximated by a quadrature rule on V. Let (g)jc denote the approximation 
of (g) , where K, is the set of 2N V + 1 indices matching the velocity grid nodes. 
By this way we compute the moments of the Maxwellian at each grid nodes 

{Xi}, 

(pi,PiUi,Ei) = {fi(t>{v))ic 

As quadrature rule we use summation over /C, providing spectral accuracy 
for smooth functions on compact support. The grid 1C is chosen to include 
most of the mass. For a given number of nodes N v , an optimal choice of 
the grid is obtained as a compromise between the extension of the velocity 
domain and the resolution of the grid. A more sophisticated strategy for the 
optimal velocity domain configuration lies in adapting the grid to the flow 
dynamics [33j.[34j. 

Once the moments are computed on the grid, they can be in turn com- 
puted in Xi, by a suitable interpolation formula, so that the Maxwellian gets 
easily evaluated. Details about the interpolation will be given in Section f3.21 



The scheme fj3.3|) can be used to perform the time step. This explicit 
scheme is first order accurate. Moreover the stability is not preserved as r 
approachs to zero (in this case the sistem (|3,3p becomes stiff). Although 
the scheme could be made higher order accurate, because of its explicit 
formulation, it requires very small time steps, becoming in turn impractical 
at the fluid regime. 

In order to circumvent the stiffness arising from the fluid regime, an 
implicit formulation of the system fj3.2j) is highly desired 



/«+ 1 = /n + A^+ 1 (/), (3.5a) 



Xi = Xi + VjAt, i = l,...,N x , j = -N v , . . . ,N V , (3.5b) 
1 

— i 

T 



H^ 1 = -(M j [fr +1 ]-f^ +1 ),. (3.5c 



Hence the moments of are needed to compute Mj[f™ +1 ], Vxj and for 

each velocity node in V. As we will see, this can be done explicitly by 
computing the moments of both sides of (|3.5ajl 

{f n+1 (x, v)<t>{v)) = (f n (x, v)tf,(v)) + (H^imv^At. (3.6) 

As pointed out in |27j , we observe that the moments of the relaxation oper- 
ator are identically zero 

(H n+ \f)<t>(v)) = ((M[P+\x,v)] - r+\x,v))ct>{v)) = 
therefore 

(f n+1 (x, v)</>(v)) = (f n (x, v)<Kv)), Vx e R. (3.7) 

and the moments can be easily computed. The Maxwellian M[f n+1 (x,v)] = 
M(v; {Pi +1 , u™ + , T™ +1 }) is known and the distribution function at the next 
time step can be explicitly evaluated 

Once Eqs. (|3.6l 13.71 13. 8p are discretized on a grid, the resultant first order 
scheme can be written as 

rh + AtM n+1 

^ ~ 7+Ai • ^ 

where M™^ = M(v; it™ +1 , T i n+1 ), and the moments are computed 
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by 



N v 

E 

j=-N v 
1 N v 

~n+T E V jfij^ V 
Pi j=-N v 

1 iV„ 

I £ #SA„-(^)7" 

This approach allows to use CFL numbers greater than one. Moreover, if 
t<1 f n (x,v) relaxes very fast toward the local Maxwellian. 
Remark Since the moments are computed by a quadrature formula, it is 
not properly true that, in the discrete formulation, M[f(t, x, v)] and f(t, x, v) 
have the same moments. To get an insight on this aspect see |23j . In that 
paper the author introduces the notion of a discrete Maxwellian, which is 
more consistent with the discrete formulation of the problem. The discrete 
BGK model obtained using such Maxwellian is conservative and entropic. By 
enough large number of grid points in velocity, the continuous and discrete 
Maxwellians give comparable results. However, for coarse discretization in 
velocity, the discrete Maxwellian introduced in |23| produces better results. 

Next section is aimed to present a more general class of methods, based 
on this basic scheme. The accuracy gets improved, preserving the discussed 
properties above. 

3.2 General WENO reconstruction 

The accuracy and the shock capturing properties of the scheme near the fluid 
regime require a suitable nonlinear reconstruction technique for the computa- 
tion of f§. ENO (Essentially-Non-Oscillstor) and WENO (Weighted-ENO) 
methods provide the desired high accuracy and non oscillatory properties 
(see |37j). Both methods are based on the reconstruction of piecewise smooth 
functions by choosing the interpolation points on the smooth side of the func- 
tion. In ENO methods these points are choosen according to the magnitude 
of the divided differences evaluated by two candidate stencils. In WENO 
methods the different polinomials, defined on the stencils, are weighted in 
such a way that the information about the function on both sides can be 
used. Here we focus on WENO reconstruction by introducing the general 
framework of the implementation. Let us consider the space grid 
and the discrete distribution function T = {fi}izz known on any space grid 
point. For simplicity the time step and velocity grid node indices are not 
used. The goal is to construct a 2n — 1 degree WENO interpolation on the 
interval Let be the Lagrange polinomial built on the stencil 



Pi 



n+1 



U 



n+1 
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n+1 
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5* — {^i— n+1 ' • ■ •) -^i+n} • 

It can be written as follow 



L(0=J2 l ™ ) (0 pik) (0> S€[xi,xn-i], k = l,...,n (3.10) 



k=l 



where lw\C) are the linear weights, of n — 1 degree polinomials, and P^ k \^) 
are n degree polinomials, interpolating J 7 on the stencil Sfc. As shown in 
Figure[2l all the stencils overlap on [x^Xj+i]. From the definition (|3.10p . the 

Si 



s„ 

Figure 2: Stencil for WENO reconstruction. 

(k) 

linear weights ly, have to vanish at the nodes outside and must have unit 
sum in the nodes of S, where they are nonzero. This condition leads us up 
to write the linear weights as follow 

4 fc) (o=7 fc n ( 3 - n ) 

x t es\s k 

where 7^ are evaluated by imposing the unit sum condition 

n 

5^4*)(xi) = l, Xi eS. (3.12) 

k=l 

For a detailed discussion about the calculation of the coefficients 7^ we refer 
to [7]. The WENO reconstruction is expressed by a simple formula 

n 

R[^](0 = J2 w ^P {k) (0- (3-13) 

k=l 

Here R [J-] (£) denotes the numerical reconstruction of the approximate so- 
lution, whereas are the nonlinear weights. To ensure stability and 
consistency the following properties are required 



Wk 



(O>0, Y j w k {x i ) = l, i,k = l,...,n. (3.14) 



k=l 



Following [37| and |21| . such conditions are satisfied by choosing 

w k (0 = k = l,...,n (3.15) 







with 

MO = T^^n (3-16) 

where the parameter e > is used in order to avoid that the denominator gets 
0, whereas are the smoothness indicators. Usually e = 10~ 6 . Following 
[37] . we use 

By this method it is possible to get high order reconstruction in specific 
points for any cell in space. 

3.2.1 Third order WENO interpolation 

In this case the stencil is S = {xj-i, Xi, Xj+i, ££+2}) an d the WENO recon- 
struction operator is given by a superposition of two parabolas pW(£) and 
P( 2 )(£), defined respectively on two overlapping stencils Si = Xj, Xj+i} 

and S 2 = {xi,x i+ i,x i+ 2}, 

R[T] (0 = w^OP^HO + MOP^HO (3-18) 
The linear weights are 

4 1) (0 = ^7 1 , 4 2) (0 = ^^, (3-19) 
and the smoothness indicators can be computed by f|3. lTj) 

_ 13 3 16 2 25 2 13 7 19 

Pi — Y2 * _1 "3" T2 ""^-./i— lji "T g Ji-lji+l "g* Jiji+l 

13 2 16 2 25 a 13 7 19 

P2 — Y^ii+2 + "^/i+l + ^2 ' 2 _ ii+2/i+l + g-Ji/i+l- 

3.3 Time discretization 



System (J372J) is a typical ordinary differential equation with relaxation, to be 
solved in the characteristics framework. Relaxation time lies in a very wide 
range. It typically extends from order one to very small values compared to 
the compared to the time scale of the problem. This is the main motivation 
leading us up to treat the relaxation operator by L-stable diagonally implict 
Runge Kutta (DIRK) schemes [THJ EH1 122] . When applied to system {32]) it 
reads 



1 1 

F®(x,v) = f n (S:,v) + AtVa.-fMpW^.)] - F^(x,v)) 

i — ' T 

k=l 

1 1 

f n+1 (x,v) = f n (x,v) + AtJ2 w k-(M[F^ k \x,v)} -F (fc) (x,v)) 



k=l 
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The triangular v x v matrix, A = (aik), and the coefficient vectors, c 



(1, ...,c„) T and w = (1, 



, are given by consistency and order conditions. 



They characterize completely a DIRK scheme, which can be rappresented by 
the Butcher's tableaux 



c 


A 




T 




W 



The internal stages are practically evaluated by a sequence of elementary 
implicit Euler steps. Scheme (|3,9p corresponds to (|3.20p with implicit Euler 
scheme. It will be denoted SI. The DIRK methods considered in this work 



arc 



S2 



a 




a 




1 


1 


— a 


a 




1 


— a 


a 



S3 



1/2 






(l + 7)/2 


(1 " 7 7)/2 


7 


1 


1-5-7 


5 7 




1-5-7 


5 7 



which are the implicit parts of the ImEx schemes ARS(2,2,2) and ARS(3,4,3) 
respectively, see |28j . and are a second and a third order L-implicit schemes. 
The coefficients are : 



a 



1 



V2 
2 ' 



7 = 0.4358665215, 5 = -0.644373171. 



The schemes constructed by the tableau above will be denoted by S2 and 
S3 respectively. The L-stability guarantees that the schemes are able to 
capture the limit r — > 0. Other choices are possible of course. The terms 
F^'(x,v) are the internal stages. We point out that the internal stages have 
to be evaluated along the characteristics solving an implicit equation, see 
Figure [4j Due to the strong nonlinearity of the relaxation operator, this 
task, in principle, is not easy to be solved. To circumvent this difficulty we 
proceed as follow. At the starting point all we have is the initial condition. 
The goal is to evaluate F± (v), in all discrete velocity domain, on the char- 
acteristics framework at the spatial coordinate x^p = X{ — v(l — c\)/S.t. The 
prepocessing calculation consists of providing a preliminary internal stage, 
which is denoted by F^(v), in each spatial grid node X{, i = 1,...,N X . 
To this end we proceed by performing a single time step of amplitude c\/\t 
using an implicit first order scheme, as shown before. Thus we can write 



and ([3 



(F^\v)4>(v)) K = {f?(v)</>(v)) K , i 1 V, 



fn {v) + ^ au[M[ F(\ v) _F(l) {v)) 



(3.21) 
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Figure 3: Characteristic field for positive velocity grid node. 



since 



((M[F^](v)-F^(vmv)) K =0. 

Once the moments, and are computed, the Maxwellian at the 

first stage M[F^](v),i = 1, ...,N X , is evaluated by 



M[fW](v 



v — u. 



(1)|2- 



exp 



f 1 l ± I 

(2wRT> ) ) 1 / 2 \ 2RT 1 



(i) 



(3.22) 



Finally the internal stage Fp~'(v) is computed by WENO reconstruction at 



the point by the values F- L, (v). 




Xi—2 Xi Xi—i 

Figure 4: Characteristic field for positive velocity grid node. 
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For multi-stage case, in order to evaluate the relaxation operator for the 
l-th internal stage 

±(M[F?\v)]-fP(v)) 

we need to compute the moments of {Fl («)}, i = 1,...,N X . To this end, 
let the set of internal stages {F^\ F^ 2 \ ...,F^~^} be already computed at 
the characteristic's set points {x^ 1 ), x^ 2 \ . . . , rc^" 1 )}, for each spatial grid 
node. Therefore we know the pre- processed stages {F^ ,F^\ ...,F ( - l ~ r >}. 
The Z-th internal stage is evaluated repeating the initial strategy, as showed 
in Figure HJ by using the pre-processed stages F^ -1 ) to get the moments 
p®,u® and and the Maxwellian M[F^ l) }{v),i = 1,...,N X . Hence it is 
possible to write 

F«\v) = f*( v ) +ail ^(M[4 l) ](v) - F«\v)) (3.23) 



where 

l-l 

ft(v) = f?(v) + At^2a lk (M[F^](v) - >(«)) 

k=l 

Equation (|3.23p is then solved for F^ l \v) by the same technique used for 
Eq. (|3.2ip . Once all internal stages are evaluated, finally we can perform the 
time evolution step 

Remark In practice the Runge Kutta fluxes can be computed from the 
internal stages. For example 

t J J an 

Hence the scheme can be used in the limit r — > 0, with no constrain on the 
time step amplitude. 



fr + Hv) = fr(v) + AtJ2w k ^(M[FP](v)-F^(v)). (3.24) 

k=l 

4 Numerical Tests 

These tests are aimed to verify the accuracy (test 1) and the shock capturing 
properties (test 2) of the schemes. 
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4.1 Test 1: regular velocity perturbation 

This test has been proposed in |27| . The solution is smooth, and the accuracy 
can be tested. Initial velocity profile is given by 

u (x) = ^ (exp (-(era; - l) 2 ) - 2exp (-(ax + 3) 2 )) , x G [-1, 1] 

where a is a positive constant parameter. Initial density and temperature 
profiles are uniform, with constant value, p = 1 and T = 1 respectively. The 
initial condition for the distribution function is the Maxwellian, computed by 
given macroscopic fields. The boundary conditions are imposed by prescribed 
moments as well. Two regimes (rarefied and fluid) have been investigated, 
corresponding to different Knundsen numbers, r = 10~ 2 and r = 10~ 6 . The 
final time, for both cases was 0.04, showed being large enough in order to get 
the thermodynamic equilibrium. Accuracy and conservation tests have been 
performed at the final time. The errors has been computed using a reference 
solution, defined on a finer grid, with N x = 1280 and N v = 20. The test case 
has been performed using N v = 20 (as for the reference solution), for each 
spatial grid nodes number, uniformly spaced in [-10,10]. Total entropy 



has been also computed by a fourth order integration formula in space (see 
|27j). for both of the Knundsen numbers and it is reported in Figure [TU1 It is 
possible to observe that the functional decreases during the time evolution, 
as expected [8] . The relative errors and order of accuracy are shown in Tables 
|1I2|3|4] . for the schemes S2 and S3. By using a reference solution at CFL = 
0.5, with N x = 200 spatial nodes (whereas N v was unchanged). Several 
computation have been carried out, for different CFL numbers. This test is 
aimed to check the correct behaviour of the schemes as the CFL changes, 
leading up to diffusity spurious errors. The results in Table \5\ show that the 
schemes work very good for each order of accuracy, by limiting the errors 
in a narrow and numerically satisfactory range. Finally the conservation 
has been investigated. Despite the schemes are not strictly conservative, 
conservation properties look good, though not many space grid nodes are 
used. That is because for smooth solutions the weights are close to their 
linear values. A conservative version of this class of schemes is introduced in 
the next section. 

4.2 Test 2: Riemann problem 

This test allows us to evaluate the capability of our class of schemes in captur- 
ing shocks, contact discontinuities and the density profile in a rarefaction. 
The macroscopic fields are initially assigned in the domain, satisfying the 
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Rankine-Hugoniot shock jump conditions. In particular we are interested to 
the behaviour in the fluid regime limit. The presented results are density, 
velocity and temperature profiles, for r = 10 -2 and r = 10 -6 , respectively. 
Like for test 1, the boundary conditions are imposed by Maxwellians com- 
puted by prescribed macroscopic moments. Total entropy calculation results 
are also presented. For this test two values r are employed , r = 10 -2 and 
r = 10 -6 . N v = 60 nodes are used in the range [-10,10] of the discrete 
velocity domain, as in |27| . 



5 The conservative version 



Here we introduce the conservative version, first order accurate. Generalizing 
to higher order of accuracy is pretty straightfoward. Let us consider the 
original problem (|2.1|) . The first internal stage F^(v) can be computed by 
solving the equation (|5.ip 



F?\v) = f?(v) - AtV(vF^) + ^(M[if >](«) - if >(«)). (5.1) 

Of course a suitable discrete form of the convective part is needed. To 
overcome this point we use a conservative finite difference approximation 
introduced in [38]. This efficient method is aimed to look for a function / 
(typically a polynomial) such that the data F^ 1 ' = vF^ are interpolated, 
in the sense of cell average 

1 r x i+i/2 A 
F = — fdx 

AX Jxi-l/2 

in such a way that 

d x F^\ Xi = ^- (/(* m/2 ) - /Wi/ 2 )) • (5.2) 

We remark that the flux depends linearly on the function F itself. The flux 
/ is splitted of course in two contributions at each cell border 

/fe+1/2) = f + ( x 7 +1/2 ) + /~( x m/ 2 )- ( 5 - 3 ) 
The reconstruction is performed on F + and F~ that are defined as follow 

f+ M". »><>, F- = {\„ " £0 ' (5.4) 
[0, v < 0, [vFW, v<0. 

Here we use a polynomial reconstruction based on WENO method, like in 
|27] (see |38| for the details). 

Let us consider the stencils {f i+ i-i, f i+ i, fi+i+i}i=-i,o,i and denote by L\{x) 
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Figure 5: Conservative diagram. 

the parabolas that interpolate the given data, in the sense of cell average. 
Since the flux is linear, the data can be identified as the distribution function 
values on the spatial nodes, as well. For each cell border we compute three 
approximations respectively for /(a^wg) an d f( x f + i/ 2 )^ as f°how 

/'( X i "+l/ 2 ) = L \( X i+l/2), f\ X t + l/ 2 ) = -^+1(^+1/2), ' = —1,0,1. 

Thus they are written as follow 

1 - 7 - 11 - 

/ _1 (^i/ 2 ) = ^ F i-2 ~ -^Fi-t—Fi 

= -1^-1 + 1% + 1**1 

1 - 5 - 1 - 
^^7+1/2) = ^Fi + -F i+ i- -F i+2 . 

The computation of f \ x 1^1/2) * s performed by using the stencils symmetri- 
cally mirrored. To compute the weights we need the smoothness indicators 

13- - -1- 

Pr = — (F i _ 2 -2F t _ 1 + F i ) 2 + -(F l _ 2 -4F i _ 1 + 3F i ) 2 , 

13 - - 1 - - 

$ = -(F l _ 1 -2F t + F t+l ) 2 + -(F l _ l -F l+1 ) 2 , 

13- - - 1-- 
Pi = -(F l -2F l+1 + F t+2 ) 2 + -(3F l -4F l+1 + F i+2 ) 2 

The nonlinear weights, w\ and w\, are computed as in (|3,15|) . but now the a 
coefficients are evaluated as follow 

a 1 = dl a 1 = j 

1 (e + Pir 1 (e + 01)* 
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where 



d-i = —, d = -, di = —, di = d_i, 1 = -1,0,1. 
10 5 10 

Finally, we have 

l i 

/>r+i/ 2 ) = E w if l ( x 7+i/2)> f (4+1/2) = E ^^«+i/ 2 )- ( 5 - 5 ) 

Summarizing, initially we perform a relaxation step, by using the first or- 
der implicit semi-Lagrangian scheme, to get a suitable predictor value of the 
internal stage, F^\v). This is used to compute f l {x^ +l j 2 ) and f l (4+1/2^ 

for each cell to get f( x ^ + i/ 2 ) an< ^ f( x i+l/2^ (|5.5p respectively. In this 
part of the computation the CFL can be greater than one, of course, since 
we know in advance the solution of the characteristic lines. Once the flux is 
evaluated at each cell border by (|5.3p . it is possible to compute the discrete 
convective term by (|5.2p . Finally the corrector step is applied by solving 
implicitly the equation (|5.ip . This procedure can be straightforwardly ex- 
tended to multiple internal stages. Figure [8] shows the comparison between 
the proposed numerical method and the reference solution computed by a 
high order solver of the Euler equations. We observe a clear improvement of 
the quality of the solution. Also in this case the accuracy, despite a slight im- 
provement, doesn't hold the third order in fluid regime, but the conservation 
is fully achieved. 



L2 — Relative errors 


N x 


Density 


Velocity 


Temperature 


20 


2.54838e-03 


2.35049e-03 


4.63423e-03 


40 


5.57339e-04 


3.64146e-04 


9.17053e-04 


80 


8.41532e-05 


5.21314e-05 


1.28062e-04 


160 


1.17817e-05 


8.94658e-06 


2.53109e-05 


320 


1.69746e-06 


1.95126e-06 


6.47027e-06 




L 


2 — Orders 




N x 


Density 


Velocity 


Temperature 


40 


2.193 


2.690 


2.337 


80 


2.727 


2.804 


2.840 


160 


2.836 


2.543 


2.339 


320 


2.795 


2.197 


1.968 



Table 1: Scheme S2, r = 10~ 2 , CFL = 4.5. 
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L2 — Relative errors 


N x 


Density 


Velocity 


Temperature 


20 


2.96809e-03 


3.15227e-03 


6.37101e-03 


40 


6.57722e-04 


6.05216e-04 


1.94043e-03 


80 


1.11120e-04 


1.36059e-04 


5.26168e-04 


160 


2.25137e-05 


4.61239e-05 


1.59816e-04 


320 


6.10643e-06 


1.63245e-05 


5.05549e-05 




L 


2 — Orders 




N x 


Density 


Velocity 


Temperature 


40 


2.174 


2.381 


1.715 


80 


2.565 


2.153 


1.883 


160 


2.303 


1.561 


1.719 


320 


1.882 


1.498 


1.660 


Table 2: Scheme S2,r = 10~ 6 , CFL = 4.5. 


L2 — Relative errors 


N x 


Density 


Velocity 


Temperature 


20 


2.41539e-03 


1.97185e-03 


4.36445e-03 


40 


4.93444e-04 


2.90747e-04 


8.14397e-04 


80 


7.36995e-05 


4.27296e-05 


1.14397e-04 


160 


1.06248e-05 


5.91413e-06 


1.54660e-05 


320 


1.55051e-06 


9.55269e-07 


2.89414e-06 




L 


2 — Orders 




N x 


Density 


Velocity 


Temperature 


40 


2.291 


2.762 


2.422 


80 


2.743 


2.766 


2.832 


160 


2.794 


2.853 


2.887 


320 


2.777 


2.630 


2.418 



Table 3: Scheme S3,r = 10~ 2 , CFL = 4.5. 

6 Conclusions 

In this work a new class of semi-Lagrangian schemes for BGK model of the 
Boltzmann equation has been introduced. These schemes show to be highly 
efficient and accurate, allowing us, not only to investigate the rarefied regime 
of the gasdynamics system, but also to converge correctly to the fluid limit 
for discontinous solution (shock wave propagation) . Moreover also high CFL 
values (less than one in classical schemes) do not affect the numerical solution 
by the diffusivity, preserving accuracy and conservation. These properties 
make the extension to multidimensional schemes highly desirable, leading 
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L2 — Relative errors 


N x 


Density 


Velocity 


Temperature 


20 


2.02024e-03 


2.72023e-03 


4.79517e-03 


40 


5.04007e-04 


7.57946e-04 


2.06197e-03 


80 


9.91878e-05 


2.63921e-04 


9.43964e-04 


160 


4.10972e-05 


1.48278e-04 


5.14779e-04 


320 


2.17256e-05 


7.52320e-05 


2.39514e-04 




L 


2 — Orders 




N x 


Density 


Velocity 


Temperature 


40 


2.003 


1.844 


1.218 


80 


2.345 


1.522 


1.127 


160 


1.271 


0.932 


0.975 


320 


1.120 


1.079 


1.104 



Table 4: Scheme S3,r = 10~ 6 , CFL = 4.5. 



CFL = 1.5 




Density 


Velocity 


Temperature 


SI 


6.000e-05 


1.120e-04 


1.700e-04 


S2 


6.000e-05 


1.120e-04 


1.700e-04 


S3 


6.000e-05 


1.120e-04 


1.700e-04 


CFL = 5.5 




Density 


Velocity 


Temperature 


SI 


2.500e-04 


2.977e-04 


4.400e-04 


S2 


2.500e-04 


2.972e-04 


4.400e-04 


S3 


2.500e-04 


2.973e-04 


4.400e-04 


CFL = 10.5 




Density 


Velocity 


Temperature 


SI 


3.000e-04 


3.120e-04 


4.700e-04 


S2 


3.000e-04 


3.115e-04 


4.700e-04 


S3 


3.000e-04 


3.117e-04 


4.700e-04 



Table 5: Comparison test for different CFL. Test 1; N x = 200, r = 10~ 2 . 

up to the possibility to simulate more realistic problems, like bidimensional 
generalized Riemann problems or microfluidics-based devices of engineering 
interest. The authors are currently go along this way, and a conservative 
multidimensional version of this class of schemes is under investigation. 
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t =le-2 


N x 


Density 


Momentum 


Energy 


20 


3.19103e-04 


6.45517e-04 


6.44489e-04 


40 


8.68364e-05 


2.43064e-05 


1.53479e-04 


80 


1.86619e-05 


6.93736e-06 


3.38482e-05 


160 


2.21369e-06 


6.03659e-07 


3.82991e-06 


320 


2.47089e-07 


6.48153e-08 


4.23732e-07 


640 


2.75503e-08 


5.92434e-09 


4.57648e-08 


1280 


3.21756e-09 


6.82768e-10 


5.29436e-09 


t =le-6 


N x 


Density 


Momentum 


Energy 


20 


3.86571e-04 


8.49545e-04 


8.59503e-04 


40 


1.25170e-04 


4.22174e-05 


2.34721e-04 


80 


3.17682e-05 


1.53056e-05 


6.50680e-05 


160 


4.65888e-06 


1.86177e-06 


9.55101e-06 


320 


5.94587e-07 


2.25448e-07 


1.20946e-06 


640 


7.47884e-08 


2.69514e-08 


1.52337e-07 


1280 


9.24536e-09 


3.27055e-09 


1.87483e-08 



Table 6: Errors in conservation. Test 1, scheme S2. 
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Figure 6: Test l,r = 10 _6 ,CFL=10.0 . Left column: from the top to the 
bottom, density, velocity and temperature for the schemes SI (cross), S2 
(circle), S3 (dashed). Right column: from the top to the bottom, density 
zoom for N x = 40, 80, 160, respectively. 
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Figure 8: Test 2. From the top to the bottom, density, velocity and pressure, 
r = 10~ 6 . Reference solution (Nx=400) (solid line) and the large time step 
DIRK scheme Nx=300, CFL=4.5 (star) 



2G 




Figure 9: Entropy for test l;scheme S3, N x = 800, r = 10 -2 (left), r = 10~ 6 
(right). 
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Figure 10: Entropy for test 2; scheme S3, N x = 800, r = 10" 2 (left), r = 
10" 6 (right). 
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